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ABSTRACT 

This  report  documents  two  Fortran  subroutines  ~  FDCALC  and  FDCORE  —  that  compute  a  set  of 
forward-difference  intervals  to  be  used  in  minimizing  a  smooth  function.  The  primary  subroutine 
FDCORE,  which  is  called  repeatedly  by  FDCALC,  produces  a  suitable  interval  for  a  univariate 
function,  as  well  as  approximations  to  the  first  and  second  derivatives. 


*FDCALC  and  FDCORE  are  available  from  the  Oflice  of  Technology  Licensing,  105  Encina  Hall, 
Stanford  University,  Stanford,  California,  04305. 

The  material  contained  in  this  report  is  based  upon  research  supported  by  the  U.S.  Department 
of  Energy  Contract  DE-AC03-768F00326,  PA  No.  DE-AT03-76ER72018;  National  Science  Foun¬ 
dation  Grants  MCS-7926000  and  KCS-8012974;  the  Oflice  or  Naval  Research  Grant  N00014-75- 
C-0287;  and  the  U.S.  Army  Research  OIRce  Contract  DAAG29-79-C-0110. 
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When  minimising  a  smooth  function  whose  derivatives  are  not  available,  it  is  common  to  use 
finite-difference  approximations  in  a  gradient- based  method.  However,  certain  “standard”  choices 
for  the  finite-difference  intervals  may  produce  unnecessarily  inaccurate  results.  The  subroutines 
documented  in  this  report  are  an  implementation  of  the  method  of  Gill  et  al.  (1983),  which  is 
designed  to  compute  “good”  intervals  for  forward-difference  approximations  to  the  gradient. 

Our  implementation  is  based  on  a  reverse-communication  control  structure.  The  “core”  sub¬ 
routine  FDCORE  (which  treats  a  single  variable)  must  be  called  repeatedly  by  an  outer  subroutine 
in  order  to  obtain  a  set  of  intervals  for  multivariate  optimisation.  The  outer  subroutine  FDGALC 
described  in  this  report  should  be  suitable  for  many  applications.  If  it  is  not,  the  user  may  develop 
an  outer  subroutine  that  caters  for  the  special  needs  of  his  problem. 

The  method  used  by  FDCORE  requires  at  least  three  evaluations  of  the  function  for  each 
component  of  the  gradient,  even  for  a  well  scaled  problem.  An  optimization  algorithm  that 
required  this  number  of  function  evaluations  at  every  iteration  would  be  uncompetitive  with 
alternative  non-derivative  methods.  Fortunately,  in  practice  several  factors  make  it  possible 
to  obtain  adequate  gradient  approximations  with  only  one  function  evaluation  per  component. 
First  and  most  important,  it  is  our  experience  that,  for  many  functions,  the  finite-difference 
intervals  generated  by  the  method  of  FDCORE  do  not  vary  significantly  from  one  iteration  to  the 
next.  Second,  these  intervals  do  not  usually  differ  widely  from  the  “optimal”  intervals.  Finally, 
finite-difference  gradient  methods  can  generally  make  satisfactory  progress  as  long  as  the  overall 
gradient  vector  has  a  reasonable  level  of  accuracy;  it  is  not  essential  for  each  component  of  the 
gradient  to  have  close- to- maximal  accuracy  at  every  iterate. 

Based  on  these  observations,  we  suggest  that  FDCALC  (or  its  equivalent)  should  be  executed 
at  a  “typical”  point  (usually,  the  initial  point  of  the  minimisation).  The  set  of  intervals  obtained 
should  then  be  used  to  compute  forward-difference  approximations  at  subsequent  iterates. 

We  emphasize  that  these  subroutines  are  not  intended  to  compute  the  most  accurate  possible 
estimate  of  the  gradient  at  a  single  point.  The  general  problem  of  finding  approximate  derivatives 
by  finite  differences  is  known  as  numerical  differentiation.  A  recent  discussion  of  methods  for 
numerical  differentiation  is  given  by  Lyncss  (1977)  (for  other  references,  sec  Gill,  Murray  and 
Wright,  1981).  Many  automatic  differentiation  routines  require  a  significant  number  of  function 
evaluations  often,  at  least  ten  per  derivative  and  hence  are  not  appropriate  within  a  finite- 
difference  gradient  method. 
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1 1.  DESCRIPTION 


We  shall  briefly  summarise  the  procedure  of  FDCORE  applied  to  a  univariate  function  /  at 
the  point  x.  (In  order  to  obtain  a  set  of  intervals  for  the  function  F(x),  where  z  is  an  n- vector, 
the  procedure  is  applied  to  each  component  of  x,  keeping  the  other  components  fixed.)  Roughly 
speaking,  the  method  is  based  on  the  fact  that  the  bound  on  the  relative  truncation  error  in 
the  forward-difference  approximation  tends  to  be  an  increasing  function  of  h,  while  the  relative 
condition  error  bound  is  generally  a  decreasing  function  of  h. 

The  "best”  interval  h,  is  given  by 


where  4  is  an  estimate  of  /"(x),  and  e A  is  an  estimate  of  a  good  bound  on  the  absolute  error 
associated  with  computing  the  function  (for  a  discussion  of  tA,  see  Chapter  8  of  Gill,  Murray  and 
Wright,  1981).  Given  an  interval  h,  4  is  defined  by  the  second-order  difference  approximation 

A  /(*  +  h)  -  2 /(»)  +  /(»  -  h) 

9  h * 

The  decision  as  to  whether  a  given  value  of  4  is  acceptable  involves  £(4),  the  following  bound 
on  the  relative  condition  error  in  4: 

(When  4  is  sero,  £(4)  is  taken  as  an  arbitrarily  large  number.) 

The  procedure  selects  the  interval  h4  (to  be  used  in  computing  4)  from  a  sequence  of  trial 
intervals  {hk}.  The  initial  trial  interval  is  taken  as  10fi,  where 

t=2,,+wVra-  (j) 

(The  quantities  rj  and  w  are  discussed  later.)  The  value  of  £(4)  for  a  trial  value  hk  is  defined  as 
“acceptable”  if  it  lies  in  the  interval  (.001,  .1].  In  this  case,  h4  is  taken  as  /»*,  and  the  current  value 
of  4  is  used  to  compute  h,  from  (1).  If  £(4)  is  unacceptable,  the  next  trial  interval  is  chosen  so 
that  the  relative  condition  error  bound  will  either  decrease  or  increase,  as  required.  If  the  bound 
on  the  relative  condition  error  is  too  large,  a  larger  interval  is  used  as  the  next  trial  value  in  an 
attempt  to  reduce  the  condition  error  bound.  On  the  other  hand,  if  the  relative  condition  error 
bound  is  too  small,  hk  is  reduced. 

The  procedure  will  fail  to  produce  an  acceptable  value  of  £(4)  in  two  situations.  Firstly,  if 
/"(*)  is  extremely  small,  then  £(4)  may  never  become  small,  even  for  a  very  large  value  of  the 
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interval.  Alternatively,  <?(♦)  may  never  exceed  .001,  even  for  a  very  small  value  of  the  Interval. 
This  usually  implies  that  /*(*)  is  extremely  large,  and  occurs  most  often  near  a  singularity. 

As  a  check  on  the  validity  of  the  estimated  first  derivative,  the  procedure  provides  a  com¬ 
parison  of  the  forward-difference  approximation  computed  with  h,  and  the  central-difference 
approximation  computed  with  h4.  If  these  values  do  not  display  some  agreement,  neither  can  be 
considered  reliable. 

The  parameters  7  and  u  in  (3)  are  local  variables  ETA  and  OMEGA  of  subroutine  FDCORE;  they 
are  both  set  to  1  with  DATA  statements,  and  may  be  changed  if  appropriate  (see  Section  10). 
Other  aspects  of  the  algorithm  may  also  be  adjusted  by  changing  the  values  of  local  variables 
in  FDCORE.  In  particular,  MAX  defines  the  maximum  number  of  trial  intervals,  BNDLO  and  BNDUP 
define  the  range  of  acceptable  values  for  C(4),  and  RHQ  defines  the  factor  by  which  hi,  is  changed 
at  each  iteration. 
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8UBR0PTINE  FDCALC (  FUN.  HSGLVL.  I.  EP8A.  X. 

NUMF,  NUMOK,  UFO. 

FX.  GRAD.  HCNTRL.  HESSD.  HFORV  > 


INTEGER  KSGLVL,  N.  NUMF.  NUMOK 

INTEGER  INFO(N) 

REAL  EPSA.  FX 

REAL  X(N) .  GRAD(N) ,  HCNTRL(N) ,  HESSD (N) .  HFORV (N) 

(The  specification  of  a  parameter  as  REAL  should  be  interpreted  as  working  precision,  which 
may  be  DOUBLE  PRECISION  in  some  circumstances.) 


•-V _v V. ••  ••  -• 
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PUI  is  a  user-provided  subroutine  that  must  define  the  function  for  which  the  gradient  is  to 
be  approximated.  FUM  must  be  declared  as  EXTERNAL  in  the  routine  that  calls  FDCALC. 
The  specification  of  FUM  is: 

SUBROUTINE  FUM  (  M.  X.  FX  ) 

INTEGER  M 

REAL  FX 

REAL  X(M) 

where 

M  is  the  number  of  variables; 

FX  should  be  set  to  the  value  of  the  function  evaluated  at  X; 

X  is  the  vector  of  N  variables  at  which  the  function  is  to  be  evaluated. 

M8GLVL  determines  the  level  of  printout  from  FDCALC.  If  KSGLVL  =  0,  there  is  no  printout 
from  FDCALC;  if  KSGLVL  =  1,  a  summary  is  printed  of  the  results  Tor  each  variable.  If 
KSGLVL  =  2,  a  full  debug  printout  is  produced  by  FDCALC  and  FDCORE. 

M  is  the  number  of  variables  (N  must  be  positive). 

EPSA  is  a  positive  number  that  should  be  a  good  bound  on  the  absolute  error  associated  with 
computing  the  function  F  at  z.  In  general,  EPSA  should  not  be  less  than  eM(l  +  |F(z)|), 
where  cM  is  the  machine  precision.  A  more  detailed  discussion  of  EPSA,  and  of  methods 
for  computing  it,  is  given  in  Chapter  8  of  Gill,  Murray  and  Wright  (1981). 

X  is  an  array  of  length  N  that  contains  the  vector  of  variables  at  which  the  set  of  intervals 

should  be  computed. 
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NUMF  gives  the  total  number  of  function  evaluations  required  to  compute  the  final  set  of 
intervals. 

NUMOK  gives  the  number  of  variables  for  which  the  intervals  were  acceptable.  If  NUMOK  =  N, 
satisfactory  intervals  were  obtained  for  all  the  variables. 

IHFO  is  an  integer  array  of  length  If  whose  j-th  component  is  set  to  the  final  value  of  INFORM 
from  FDCORE  for  the  j-th  variable.  The  meaning  of  each  value  is  explained  under  INFORM 
in  Section  8. 

IX  is  the  value  of  the  function  evaluated  at  the  input  vector  X. 

GRAD  is  an  array  of  length  N  whose  j-th  component  contains  the  best  estimate  of  the  first 

partial  derivative  for  the  j-th  variable. 

HCMTRL  is  an  array  of  length  N  whose  j-th  component  is  the  best  interval  found  for  computing 
a  central-difference  approximation  to  the  partial  derivative  for  the  j-th  variable. 

BB88D  is  an  array  of  length  N  whose  j-th  component  is  the  estimate  of  the  second  partial 
derivative  with  respect  to  the  j-th  variable  (i.e.,  the  j-th  diagonal  element  of  the  Hessian 
matrix  V2F(z)). 

HFQRV  is  an  array  of  length  N  whose  j-th  component  is  the  best  interval  found  for  computing 
a  forward-difference  approximation  to  the  partial  derivative  for  the  j-tn  variable. 


SPECIFICATION  OF  FU.ORE 


FDCORE/T 


6 *  SPECIFICATION  OF  FDCORE 

SUBROUTINE  FDCORE (  DEBUG,  EPSA,  FX.  X. 

IENTRY,  INFORM.  ITER, 

CDEST,  CD SAVE,  ERRBND,  FBACK,  FDOPT, 

FDSAVE,  FFORW,  H.  HOPT,  HPHI,  HSAVE, 

OLDH.  SDEST,  SDSAVE  ) 

DEBUG 

IENTRY.  INFORM.  ITER 
EPSA.  FX,  X. 

CDEST.  CDSAVE.  ERRBND.  FBACK.  FDOPT, 

FDSAVE,  FFORW,  H,  HOPT.  HPHI.  HSAVE. 

OLDH.  SDEST.  SDSAVE 

(The  specification  of  a  parameter  as  REAL  should  be  interpreted  as  working  precision,  which 
may  be  DOUBLE  PRECISION  in  some  circumstances.) 


LOGICAL 

INTEGER 

REAL 


*  This  section  need  not  be  read  by  those  calling  FDCALC  directly. 
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DEBUG  is  a  logical  variable  that  should  be  set  to  .TRUE,  if  a  detailed  printout  is  desired.  If 
DEBUG  is  .FALSE.,  there  is  no  printout  from  FDCORE. 

EPSA  is  a  positive  number  that  should  be  a  good  bound  on  the  absolute  error  associated  with 
computing  the  function  at  X.  In  general,  EPSA  should  not  be  less  than  eM(l  4-  |/(x)|), 
where  eM  is  the  machine  precision.  A  more  detailed  description  of  EPSA,  and  techniques 
for  computing  it,  are  given  in  Chapter  8  of  Gill,  Murray  and  Wright  (1981).  The  value 
of  EPSA  must  not  be  changed  until  after  the  final  exit  from  FDCORE  for  a  given  point. 

FX  must  be  set  to  the  value  of  the  function  at  the  point  where  the  derivative  is  to  be 

approximated.  The  value  of  FX  must  not  be  altered  until  after  the  final  exit  from 
FDCORE  for  a  given  point. 

X  must  be  set  to  the  value  of  x  (the  point  where  the  derivative  is  to  be  approximated). 

The  value  of  X  must  not  be  altered  until  after  the  final  exit  from  FDCORE  for  a  given 
variable. 


*  This  section  need  not  be  read  by  those  calling  FDCALC  directly. 
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Many  of  the  following  parameters  are  included  in  the  calling  sequence  of  FDCORE  simply  to 
allow  communication  between  iterations.  For  such  a  parameter,  the  description  below  includes 
only  the  meaning  of  the  parameter  on  exit  from  FDCORE. 

IENTRY  controls  the  logic  of  FDCORE.  On  entry  to  FDCORE,  IENTRY  designates  the  portion  of 
FDCORE  to  be  executed.  IENTRY  must  be  set  to  zero  before  the  first  call  of  FDCORE  for 
each  variable.  Thereafter,  IENTRY  should  not  be  altered  by  the  calling  routine.  On 
exit  from  FDCORE,  IENTRY  indicates  the  action  to  be  taken  by  the  calling  routine,  and 
thus  the  value  of  IENTRY  should  be  tested  after  each  call  of  FDCORE.  If  IENTRY  =  5  on 
exit  from  FDCORE,  the  procedure  has  terminated  (and  the  parameter  INFORM  described 
below  indicates  the  reason  for  termination).  If  IENTRY  <  4,  FDCORE  should  be  called 
again.  If  IENTRY  =  4,  a  new  value  must  be  assigned  to  FFORV  before  calling  FDCORE.  If 
IENTRY  <  4,  new  values  must  be  assigned  to  FFORV  and  FBACK  before  calling  FDCORE. 

INFORM  indicates  the  final  result  of  FDCORE,  as  follows. 

Value  Definition 

0  The  algorithm  terminated  successfully.  The  forward-difference  es¬ 

timate  FDOPT  and  the  central-difference  estimate  computed  with  HPHI 
agree  to  at  least  half  a  decimal  digit. 

1  The  function  appears  to  be  constant.  The  variable  HOPT  is  set  to  the 
value  h  (equation  (3))  corresponding  to  a  well  scaled  problem,  and  HPHI 
is  set  to  10/»;  FDOPT,  SDEST  and  ERRBND  are  set  to  tero. 

2  The  function  appears  to  be  linear  or  odd.  The  variables  HOPT  and  HPHI 
are  set  to  the  smallest  interval  with  acceptable  bounds  on  the  rela¬ 
tive  condition  error  in  the  forward-  and  backward-difference  estimates. 
The  variable  FDOPT  is  set  to  the  corresponding  forward-difference  es¬ 
timate,  and  SDEST  is  set  to  tero. 

3  The  second  derivative  of  the  function  appears  to  be  so  large  that 
it  cannot  be  reliably  estimated.  The  variables  HOPT  and  HPHI  are 
set  to  the  smallest  trial  interval;  FDOPT  and  SDEST  arc  set  to  the 
corresponding  finite-difference  estimates. 

*  This  section  need  not  be  read  by  those  calling  FDCALC  directly. 
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The  algorithm  terminated  with  an  apparently  acceptable  estimate  of 
the  second  derivative.  However,  the  forward-difference  estimate  with 
HOPT  and  the  central-difference  estimate  with  HPHI  do  not  agree  to  half 
a  decimal  place.  This  value  of  INFORM  usually  occurs  when  the  first 
derivative  is  small.  In  this  case,  although  the  derivative  estimates  will 
probably  have  poor  relative  accuracy,  the  interval  is  still  likely  to  be 
acceptable. 

ITER  is  the  count  of  the  iteration  within  FDCORE. 

GDEST  is  usually  the  value  of  the  central-difference  derivative  estimate  with  the  input  value  of 
H.  When  IENTRY  =  4  on  entry  to  FDCORE,  the  value  of  CDEST  is  not  altered. 

GDSAVE  is  the  central-difference  estimate  computed  with  HSAVE  if  HSAVE  is  positive.  Otherwise, 
CDSAVE  is  undefined. 

ERRBND  is  a  bound  on  the  estimated  error  in  the  final  forward-difference  approximation.  (When 
INFORM  =  1,  ERRBND  is  set  to  sero.) 

FBACK  must  be  defined  by  the  calling  routine  as  follows.  If  IENTRY  <  4  on  exit  from  FDCORE, 
FBACK  must  be  set  to  the  value  of  /(X  —  H)  before  the  next  call  of  FDCORE.  FBACK  need 
not  be  defined  by  the  calling  routine  when  IENTRY  =  0  on  entry  to  FDCORE,  or  when 
IENTRY  —  4  or  5  on  exit  from  FDCORE. 

FDOPT  is  the  “best”  forward-difference  estimate  of  the  first  derivative  after  the  final  exit  from 
FDCORE.  If  INFORM  =  0  or  4,  FDOPT  is  the  forward-difference  estimate  computed  with 
HOPT.  Otherwise,  FDOPT  is  defined  as  described  above  under  INFORM. 

FDSAVE  is  the  value  of  the  forward-difference  approximation  computed  with  HSAVE  when  HSAVE 
is  positive.  Otherwise,  FDSAVE  is  undefined. 

FFORV  must  be  defined  by  the  calling  routine  as  follows.  If  IENTRY  <  5  on  exit  from  FDCORE, 
FFORV  must  be  set  to  the  value  of  /(X  +  H)  before  the  next  call  of  FDCORE.  FFORV  need 
not  be  defined  by  the  calling  routine  when  IENTRY  =  0  on  entry  to  FDCORE,  or  when 
IENTRY  =  5  on  exit  from  FDCORE. 

H  is  the  current  value  of  the  difference  interval.  On  entry  to  FDCORE  (after  the  first  call), 

H  is  the  interval  for  which  the  new  function  values  have  been  evaluated.  On  exit  from 
FDCORE,  H  is  the  interval  for  which  new  function  values  should  be  computed. 

HOPT  is  the  final  estimate  of  the  “optimal”  forward-difference  interval  when  FDCORE  terminates 
with  INFORM  =  0  or  4.  The  values  assigned  to  HOPT  for  other  values  of  INFORM  are 
described  above  under  INFORM. 
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HPHI 

REAVE 

OLDH 

SDEST 

SDSAVE 


is  the  interval  used  to  compute  the  final  estimate  of  the  second  derivative  when  FDC08E 
terminates  with  INFORM  =  0  or  4.  The  values  assigned  to  HPHI  for  other  values  of 
INFORM  are  described  above  under  INFORM. 

is  set  to  —1  in  FDCORE  when  I ENTRY  =  0.  Thereafter,  a  positive  value  of  B8AVE  on 
exit  from  FDCORE  is  the  smallest  interval  for  which  the  bounds  on  the  relative  condition 
error  in  the  forward-  and  backward-difference  approximations  are  acceptable. 

is  the  value  of  H  from  the  previous  call  of  FDCORE. 

is  the  best  estimate  of  the  second  derivative  (the  one  used  to  compute  HOPT)  when 
INFORM  =  0  or  4.  The  values  assigned  to  SDEST  for  other  values  of  INFORM  are  described 
above  under  INFORM. 

is  the  estimate  of  the  second  derivative  computed  with  HSAVE  if  HSAVE  is  positive  on 
exit  from  FDCORE. 
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10.  ERROR  RECOVERY 


Whenever  the  value  of  INFORM  from  FDCORE  is  non-sero  for  no  apparent  reason,  the  user  should 
re-run  FDCORE  with  the  parameter  DEBUG  set  to  .TRUE.,  in  order  to  obtain  a  detailed  printout  of 
each  iteration  within  FDCORE.  The  main  quantities  of  interest  are  FDCERR,  the  estimated  condition 
error  in  the  forward-difference  approximation,  and  SDCERR,  the  estimated  condition  error  (2)  in 
the  second-derivative  estimate. 


Termination 
INFORM  =  1 


INFORM  =  2 


INFORM  =  3 


INFORM  “  4 


Recommended  Action 

This  value  occurs  when  the  estimated  relative  condition  error  in  the  first 
derivative  approximation  is  unacceptably  large  for  every  value  of  the 
finite-difference  interval.  If  this  happens  when  the  function  is  not  con¬ 
stant,  the  initial  interval  may  be  too  small;  in  this  case,  the  constants  ETA 
and  OMEGA  in  (3)  may  be  changed  in  the  DATA  statements  of  FDCORE.  This 
error  may  also  occur  if  the  function  evaluation  includes  an  inordinately 
large  constant  term,  or  if  EPSA  is  too  large. 

In  this  case,  the  estimated  relative  condition  error  in  the  second  deriva¬ 
tive  approximation  remained  large  for  every  trial  interval,  but  the  es¬ 
timated  error  in  the  first  derivative  approximation  was  acceptable  for  at 
least  one  interval.  This  usually  means  that  the  function  is  linear  or  odd. 
If  it  is  not,  the  user  should  check  whether  the  estimated  relative  condi¬ 
tion  error  in  the  second  derivative  (SDCERR  in  the  printout)  appears  to  be 
decreasing  as  the  trial  intervals  increase.  If  so,  it  may  be  worthwhile  to 
alter  the  ETA  and  OMEGA  in  (3)  so  that  the  initial  trial  interval  is  larger, 
or  to  allow  more  iterations  by  increasing  UMAX  (see  Section  2). 

This  value  occurs  when  the  relative  condition  error  estimate  in  the 
second  derivative  remained  very  small  for  every  trial  interval.  This 
usually  occurs  when  the  second  derivative  is  extremely  large  (e.g.,  near 
a  singularity).  If  this  is  not  the  case,  the  user  should  check  whether  the 
value  of  SDCERR  in  the  debug  printout  is  increasing  as  the  trial  intervals 
decrease.  If  so,  it  may  be  worthwhile  to  alter  ETA  and  OMEGA  so  that  the 
initial  interval  is  smaller,  or  to  increase  KMAX  (see  Section  2).  This  error 
may  also  occur  when  the  given  value  of  EPSA  is  not  a  good  estimate  of 
a  bound  on  the  absolute  error  in  the  function  (i.e.,  EPSA  is  too  small). 

The  usual  reason  that  the  forward-  and  central-difference  estimates  fail 
to  agree  is  that  the  first  derivative  is  small.  If  the  first  derivative  is  not 
small,  it  may  be  helpful  to  execute  the  procedure  at  a  different  point. 
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II.  IMPLEMENTATION  INFORMATION 


Both  FDCORE  and  FDCALC  have  been  written  in  ANSI  (1966)  Fortran  and  tested  on  an  IBM 
3081  computer  using  the  WATFIV  Compiler  (Version  1,  Level  6).  All  subroutines  are  PFORT- 
compatible  (Ryder,  1974).  Before  the  first  call  of  FDCALC  or  FDCORE,  the  subprogram  MCHPAR 
must  be  called  to  assign  various  machine-dependent  parameters.  These  parameters  are  stored  in 
the  array  VMACH(15),  which  is  stored  in  the  labelled  COMMON  block  SOLMCH. 

The  version  of  MCHPAR  provided  by  the  Systems  Optimisation  Laboratory  contains  the  parar 
meters  associated  with  double  precision  on  a  machine  in  the  IBM  370  series  (a  listing  of  the  IBM 
version  of  MCHPAR  is  given  in  Section  12).  The  user  must  substitute  a  version  of  MCHPAR  that  is 
sppropriste  for  the  machine  to  be  used.  The  specification  of  MCHPAR  is 


SUBROUTINE  MCHPAR 

REAL  VMACH 

COMMON  /SOLMCH/  VMACH (16) 


The  first  eleven  components  of  the  REAL  array  VMACH  must  be  set  in  MCHPAR.  The  components 
of  VMACH  are  defined  as  follows. 

Value  Definition 


VMACH(l) 

VMACH(2) 

VMACH(3) 

VMACH(4) 

VMACH(5) 

VMACH(8) 

VMACH(7) 

VMACH(8) 

VMACH(9) 

VMACH(IO) 

VMACH(ll) 


is  NBASE,  the  base  of  floating-point  arithmetic. 

is  MDIGIT,  the  number  of  NBASE  digits  of  precision. 

is  EPSMCH,  the  floating-point  precision. 

is  RTEPS,  the  square  root  of  EPSMCH. 

is  FLMIN,  the  smallest  positive  floating-point  number. 

is  RTMIN,  the  square  root  of  FLMIN. 

is  FLMAX,  the  largest  positive  floating-point  number. 

is  RTMAX,  the  square  root  of  FLMAX. 

is  UNDFLV,  which  specifies  whether  or  not  underflow  is  checked  for  in 
certain  computations  (not  relevant  to  FDCORE). 

is  MIN,  the  file  number  for  the  input  stream. 

is  NOUT,  the  file  number  for  the  output  stream. 
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The  values  of  XBASE,  MDIGIT,  EPSMCH,  FLMIM  and  FLMAX  for  several  machines  are  given  in  the 
following  table,  for  both  angle  and  double  precision;  RTEPS,  RTMIM  and  RTMAX  may  be  computed 
using  Fortran  statements.  The  values  XIX  and  XQUT  depend  on  the  machine  installation. 

For  each  precision,  we  give  two  values  for  EPSMCH,  FLMIM  and  FLMAX.  The  first  value  is  a 
Fortran  decimal  approximation  of  the  exact  quantity;  use  of  this  value  in  MCHPAR  should  cause 
no  difficulty  except  in  extreme  circumstances.  The  second  value  is  the  exact  mathematical 
representation. 


Table  of  machine-dependent  parameters 


Variable 

IBM  360/370 
Single 

CDC  6000/7000 
Single 

DEC  10/20 
Single 

Univac  1100 
Single 

DEC  VAX 
Single 

XBASE 

16 

2 

2 

2 

2 

MDIGIT 

6 

48 

27 

27 

24 

EPSMCH 

9.54E-7 

16"* 

7.11E-15 

2-« 

7.46E-9 

2"*t 

1.50E-8 

2-*s 

1.20E-7 

2“** 

FLMIX 

_ 

1. OB-78 

16"*» 

1.0E-293 

2-S7S 

1.0E-38 

2-l»9 

1.0E-38 

2-ia» 

1.0E-38 

2-1*8 

FLMAX  '  1.0E+75 

:  I6M(1  !«-•) 

1. OE+322 
21070(1  _2~4«) 

1 .08+38 
2l*7(l.2-»7} 

1.0E+38 

2187(1_2-87) 

1.0E*38 

2187(l_2-*4) 

Variable 

IBM  360/370 
Double 

CDC  6000/7000 
Double 

DEC  10/20 
Double 

Univac  1100 

Double 

DEC  VAX 

Double 

XBASE 

16 

2 

2 

2 

2 

MDIGIT 

14 

96 

62 

61 

56 

EPSMCH 

2.22D-16 

16"u 

2.53D-29 

2~ 95 

2 . 17D-19 

2~99 

8.68D-19 

2-eo 

2.78D-17 

2-55 

FLMIM 

1.0D-78 

16-** 

1.0D-293 

2-*T9 

1.0D-38 

2 -179 

1.0D-308 

2~ioas 

1.0D-38 

2-198 

FLMAX 

1 .OD+75 
16m(1-18"h) 

1 .OD+322 
2,070(i_2-99j 

1 . OD+38 
2187(1_2-«9) 

1.0D+307 

2I08S(1_2-61) 

1 .00*38 

2m(l-2"®#) 
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IS.  EXAMPLE  PROGRAM  AND  OUTPUT 


This  section  contains  a  listing  and  the  computed  results  from  a  sample  main  program  that 
calls  PDCALC  to  obtain  a  set  of  intervals  for  the  following  four- variable  function 

F{*)  =  2 

»— 1 


where 

/i(l)  =  21* +  41; 

/a(0  =  exp(lOi); 

/,(<)  =  10-4!f  +  l; 

/*(<)  =  21*  -  2.51s  - 1. 

(Note  that  F  is  the  sum  of  univariate  functions.)  The  analytic  first  and  second  derivatives  of  the 
functions  {/,•}  are: 

/{  =  81s  +  4,  /f-  wi; 

fi  =  10  exp  (101),  fS  =  100  exp  (101); 

fi  =  .0002t  +  1,  f!  —  .0002; 

/'  =  81s  -  51  -  1,  /?=  121-5. 

The  point  at  which  derivatives  are  to  be  estimated  isi=(l,  .25,  10,  1  +  y/tZ)T,  where 
is  the  machine  precision.  The  corresponding  values  of  F  and  {/«}  (to  fifteen  figures)  are 

F{x)  =  26.6924930607035; 

/i(*0=  6.0;  /8(xa)  =  12.1824039607035; 

/,(*,)=  10.01;  /4(x4)  =  -1.5. 

The  exact  gradient  at  this  point  (to  15  figures)  is 

VF(*)  =  (  10,  121.82403960703,  1.002,  7y/tZ  +  <k*)r, 

and  the  exact  Hessian  diagonals  are 

12,  1218.2403060703,  2  X  10“4,  and  7  + 12^. 

The  computed  results  were  obtained  on  an  IBM  3081  computer  with  double  precision  arith¬ 
metic  (eM  »  2.2  X  10“,#). 
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1  C  whhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhmhhw 

2  C  SAMPLE  MAIN  PROGRAM  FOR  TESTING  SUBROUTINES  FDCALC  AND  FDCORE 

3  C  WITH  4-VARIABLE  EXAMPLE. 

4  C 

5  C  DOUBLE  PRECISION  VERSION  t.1.  JUNE  IMS. 

4  C 

7  C  SYSTEMS  OPTIMIZATION  LABORATORY,  DEPARTMENT  OF  OPERATIONS  RESEARCH 

6  C  STANFORD  UNIVERSITY,  STANFORD,  CALIFORNIA  94305. 

9  C 

to  C  <hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhw 
ft  c 

12  INTEGER  HSGLVL,  N,  NUMF,  NUHOK 

13  INTEGER  INFO! 4) 

14  DOUBLE  PRECISION  EPSA,  EPSHCH,  FX 

15  DOUBLE  PRECISION  6RAD14),  HCNTRLI4I,  HESSD14),  HF0RN14),  X14) 
14  C 

17  C  THE  ARRAY  NMACH1 15)  IN  THE  LABELLED  COtMN  AREA  SOLMCN  CONTAINS 
IB  C  MACHINE-DEPENDENT  RUANTITIES  NEEDED  BY  FDCALC  AND  FDCORE. 

19  C  THE  HMACH  ARRAY  IS  INITIALIZED  BY  CALLING  THE  SUBROUTINE  HCHPAR. 

20  C 

21  DOUBLE  PRECISION  NMACH 

22  COMMON  /SOLMCH/  WKACHt 15) 

23  EXTERNAL  FUN 

24  C 

25  DOUBLE  PRECISION  OABS,  DSRRT 

24  C 

27  C  THE  SUBROUTINE  MCHPAR  INITIALIZES  THE  M1ACN  ARRAY  IN  THE 

28  C  LABELLED  COMMON  AREA  SOLMCH. 

29  C 

30  CALL  MCHPAR 

31  EPSMCH  =  NMACH  ( 3 ) 

32  C 

33  N  =4 

34  C 

35  C  HSGLVL  *  1  HILL  GIVE  A  SUMMARY  PRINTOUT  FOR  EACH  VARIABLE. 

34  C 

37  HSGLVL  =  1 

38  C 

39  C  SET  THE  POINT  AT  WIICH  THE  DERIVATIVES  ARE  TO  BE  ESTIMATED. 

40  C 

41  X(1)  s  1.00*0 

42  X(2)  s  .250*0 

43  X(3)  «  10.00*0 

44  X14)  3  1.00*0  *  0SRRT1 EPSMCH ) 

45  C 

44  C  COMPUTE  EPSA,  A  BOUND  ON  THE  ABSOLUTE  PRECISION  OF  THE  FUNCTION. 

47  C 

48  CALL  FUN  1  N,  X,  FX  ) 

49  EPSA  s  10. 0D*0*EPSMCH*( 1 . 00*0  *  OABS(FX)) 

50  C 

51  CALL  F0CALC1  FUN,  HSGLVL,  N,  EPSA,  X, 

52  «  NUMF,  NUMOK,  INFO, 

53  «  FX,  GRAD,  HCMTRL,  HES30,  HFORM  ) 

54  C 

55  STOP 
54  C 

57  C 

58  C  END  OP  MAIN. 

59  END 
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605  SUBROUTINE  MCHPAR 

606  C 

607  DOUBLE  WtECIOION  WUN 

600  COMMON  /SOUCH/  MMACMI 15) 

600  C 

610  C  MCHPAR  HUOT  DEFINE  THE  RELEVANT  MACHINE  PARAMETER*  AS  FOLLflM. 

611  C  NNACH(I)  ■  NBASE  >  BASE  OP  FLOATINO- POINT  ARITHMETIC. 

612  C  NMACHI2)  *  NDI6IT  a  NO.  OF  BASE  W1ACM( 1 )  DIBITS  OP  PRECISION. 

613  C  I0UCHC3)  >  EPSMCH  ■  FLOATINB” POINT  PRECISION. 

614  C  NMACH14I  «  STEPS  «  SORT(EPSMCH). 

615  C  M1ACH15)  ■  FLMIN  «  SMALLEST  POSITIVE  FLOATINB- POINT  NUMBER. 

616  C  M1ACH16)  ■  RTMIN  a  SORT  (FLMIN). 

617  C  HMACH17)  *  FLMAX  «  LARGEST  POSITIVE  FLOATINB- POINT  NUMBER. 

618  C  I01ACHC8)  >  RTMAX  «  SORT) FLMAX). 

619  C  NMACH(O)  >  UNDFLH  *  0.0  IF  UNDERFLOM  IS  NOT  FATAL.  *VE  01HERHISE. 

620  C  NMACHMO)  «  NIN  a  STANDARD  FILE  NUMBER  OP  THE  INPUT  STREAM. 

621  C  W1ACHM1)  a  NOUT  a  STANDARD  PILE  NUMBER  OP  THE  OUTPUT  STREAM. 

622  C 

623  INTEGER  NBASE.  NDIBZT.  MM.  NOUT 

624  DOUBLE  PRECISION  OSORT 

625  C 

626  NBASE  «  16 

627  NDIGIT  a  14 

620  HMACNC 1 )  *  NBASE 

629  WIACHIC)  a  NDIBZT 

630  MIACM(3)  *  MMACH(1)<H»(1  -  NDI6IT) 

631  HMACH14)  «  DSGRTCNMACHt 3)) 

632  W1ACH( 5 )  *  MMACH< 1 )**(-62l 

633  HMACHI6)  a  DSORTI HNACM1 5 ) ) 

634  NMACN17)  a  NMACH(U**61 

635  WIACH(S)  a  DSGRT(NhACH(  7) ) 

636  M1ACM19)  a  0.00*0 

637  NIN  aff 

635  NOUT  a  6 

639  NMACHMO)  «  NIN 

640  NMACMM1  )  a  NOUT 

641  C 

642  C - M  NATPIV.  ALLOW  UP  TO  100  UNDERFLOWS. 

643  c - CALL  TRAPS  (  0.0.100  ) 

644  RETURN 

645  C 

646  C  END  OP  HCHPAR 

647  END 
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*****  OUTPUT  FROM  FDCALC.  DETAILS  FOR  VARIABLE  1 

INFORM  AND  ERROR  BOOM)  0  1.7179940-06 

ESTIMATED  GRADIENT  AND  HESSIAN  DIAGONALS  1.0000000  01  1.2000000  01 

BEST  INTERVALS  FOR  FORWARD  AND  CENTRAL  DIFFERENCES  1 .4316620-07  1.8848640-06 

FDCALC  NEEDED  3  FUNCTION  EVALUATIONS  FOR  THIS  VARIABLE 

*****  OUTPUT  FROM  FDCALC.  DETAILS  FOR  VARIABLE  2 

INFORM  AMI  ERROR  BOUN)  0  1. 7310470-05 

ESTIMATED  GRADIENT  AND  HESSIAN  DIAGONALS  1 .2182490  02  1.216304D  03 

BEST  INTERVALS  FOR  FORWARD  AND  CENTRAL  DIFFERENCES  1. 42086 70-08  1.1780400-07 

FDCALC  NEEDED  S  FUNCTION  EVALUATIONS  FOR  THIS  VARIABLE 

*****  OUTPUT  FROM  FDCALC.  DETAILS  FOR  VARIABLE  3 

INFORM  AND  R1R0R  BOUND  0  7.0136830-09 

ESTIMATED  GRADIENT  AND  HESSIAN  DIAGONALS  1.002000000  2.0000000-04 

BEST  INTERVALS  FOR  FORWARD  AND  CENTRAL  DIFFERENCES  3.5068420-05  1.0366750-03 

FDCALC  NEEDED  7  FUNCTION  EVALUATIONS  FOR  THIS  VARIABLE 

*****  OUTPUT  FROM  FDCALC.  DETAILS  FOR  VARIABLE  4 

INFORM  AND  ERROR  BOUND  4  1.3120460-06 

ESTIMATED  GRADIENT  AND  HESSIAN  DIAGONALS  7.580661D-07  6.999000D  00 

BEST  INTERVALS  FOR  FORWARD  AND  CENTRAL  DIFFERENCES  1 .8746200-07  1 .8848640-06 

FDCALC  NEEDED  3  FUNCTION  EVALUATIONS  FOR  THIS  VARIABLE 
ABNORMAL  EXIT  FOR  THIS  VARIABLE  BECAUSE.... 

THE  FORWARD  AND  CENTRAL  ESTIMATES  ARE  NOT  CLOSE. 


TOTAL  NUMER  OF  FUNCTION  EVALUATIONS  ■ 


19 


NUMOK  « 


3 
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